Finite element analysis of stresses on unsymmetric quadrilateral meshes: numerical errors and accuracy evaluation.
Since the finite element analysis is an approximate method, the question of its accuracy has always been a critical issue. On each particular problem the most obvious procedure would be to perform several FEA runs on one and the same model meshed with gradually decreasing elements’ sizes until a required accuracy level is achieved.

We decided to investigate the FEA accuracy on a particular example for which an exact analytical solution is known. The model, however, has some additional features that can provide an interesting alternative basis for the accuracy evaluation.
Our ultimate goal will be accurate assessment of the stress values. The FEA is based on nodal displacements, and stresses are essentially derivatives of the displacements (subject to the material constants) so their modeling accuracy is always questionable. If, for example, standard 4-noded elements comprise a mesh, displacements over each particular element are bilinear functions. The stress tensor components, therefore, are linear function in single direction. Since the actual stress solutions must be continuous functions over the whole model, interpolation techniques must be used, and the accuracy assessment problem becomes very complex.

Consider a rectangular plate with small circular hole in the center. In a theoretical case when the plate has infinite width the analytical solution is well known.

According to Timoshenko and Goodier (S.P.Timoshenko, J.N.Goodier: “Theory of Elasticity”, 3rd ed. McGraw-Hill Book Company), if a distributed constant orthogonal load S is applied to the plate’s left and right borders, the stress tensor components in polar coordinates are:
The most important (maximum) stresses are the tensile components in points A and B, i.e at -90 and +90 degrees correspondingly: the analytical solution gives us the normal stress value of 3S.
Consider a finite element model close to the above. We set the hole diameter significantly smaller than the plate actual dimensions, so an influence of the plate’s limited width would be rather small.

The first important property of our model would be that we use 9-noded Lagrangian elements; this guarantees high numerical accuracy. The plate will have the following parameters:
We will conduct several numerical experiments on meshes with different geometrical properties.

  1. Let the mesh be completely symmetric with respect to the vertical axis, as well as with respect to the horizontal axis, with the hole diameter equal 6 mm as before. Set the step size = 2 mm on all sides and along the two division lines along OY axis.

    The solution is expectedly symmetrical, ie. the corresponding stresses in points A and B are equal to each other. The actual difference between numerical and analytical solutions is about 4%; this is represented by diagram and table below.








  1. Now we increase the vertical step size below the hole to 2.4444444… , ie. the number of steps on the bottom part of the axis of symmetry is reduced to 18.

    The stress results are provided in the table below (“+” correspond to the control points above the OX axis; “-“ to the points below OX):


    With respect to the horizontal axis OX the meshes are very different; nevertheless, the obtained solutions are almost coincident. Here we need to emphasize that the numerical discrepancy is primarily due to the fact that the stress gradients become very large in the vicinity of the hole, so the mesh with larger steps delivers less accurate results. From practical point of view, the differences in stress values are negligible.

  2. Finally we conduct the most important experiment: we will use a mesh deliberately created with a significant degree of irregularity. That mesh will have three major properties:

    • the elements’ sizes will be “approximately” the same as in the above two experiments;
    • in the most important area of the model, ie. around the hole, the mesh has different step sizes: there are 8 steps on the right and 10 steps on the left;
    • the mesh is not symmetric with respect to any axis.




One could have expected that the numerical solution would show significant variations from each other in the points where they must be equal. The reality is the opposite:
Regardless of whether the mesh is symmetric or not, the solutions are practically equal: the difference between stresses in points A and B is about 0.5%.



This result has two very significant practical consequences:

  1. In practical FEA modeling, there is no need whatsoever to impose a requirement that on geometrically symmetrical objects symmetric meshes only must be generated. The only requirement is mesh quality, i.e. whether the mesh doesn’t contain ill-conditioned elements. That is usually achieved by selecting step size(s) along the boundary sufficiently small, to give the algorithm a greater “degree of freedom”.
  2. If on one and the same model with symmetric loads an unsymmetric mesh is generated, just a single run of the FEA modeling program is sufficient to obtain a reliable estimate of the actual numerical error (since the FEA is an approximate method by itself). This gives an alternative approach to running a FEA modeling program many times with gradually decreasing element sizes with the purpose to estimate the numerical solution quality. The suggested approach is especially valuable because for very large models of millions elements each run may require significant time, whereas often the results must be obtained as soon as possible.
All meshes used in the numerical experiments above have been generated by quadrilateral mesh generation software developed by Computational Mechanics Australia Pty.Ltd.

Copyright(c) Computational Mechanics Australia Pty Ltd - All Rights Reserved ABN 39 081 999 135